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The  finite  element  method  is  applied  to  several  simple  cases  of 
steady  flow  of  a  perfect,  incompressible  fluid.  It  is  shown  that  the 
finite  element  representation  accurately  reflects  the  behavior  of  the 
classical  flow  equations.  Some  finite  element  theory  is  included  for 
the  reader  who  may  be  unfamiliar  with  this  subject. 
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SECTION  I 
INTRODUCTION 


Finite  elements  form  the  basis  for  a  versatile  analysis  procedure  applicable  to  problems 
in  several  different  fields.  Earliest  applications  were  to  problems  in  structural  mechanics 
(Reference  1).  In  recent  years  non- structural  problems  have  also  been  treated  by  this 
method  (Reference  2). 


Basically,  the  finite  element  method  represents  an  approximate  procedure  for  satisfying 
the  problem  in  terms  of  its  variational  formulation.  In  structural  mechanics,  for  example, 
this  is  generally  accomplished  by  determining  displacement  fields  based  on  satisfying  the 
minimum  potential  energy  theorem.  Consequently,  finite  elements  furnish  a  useful  alternative 
scheme  for  applying  the  well-known  Ritz  method.  In  developing  the  method  for  nonstructural 
problems,  it  is  essential  that  the  appropriate  variational  expressions  be  known  beforehand. 
For  the  flow  problems  taken  up  in  this  paper,  such  expressions  are  well  known. 

The  area  or  volume  defining  the  total  region  of  a  problem  may  be  considered  to  be  de¬ 
composed  into  an  assemblage  of  subregions.  For  example,  in  the  case  of  a  two  dimensional 
flow  these  subregions  may  be  taken  as  triangles  of  arbitrary  size  and  shape.  Each  such 
triangular  region  represents  a  finite  element.  Vertices  of  the  triangles  become  nodal  points 
of  the  finite  element  assemblage.  They  assume  a  particular  importance  in  the  development 
of  the  finite  element  equations.  Figures  1  and  2  show  typical  assemblages  based  on  the 
triangular  element. 

The  governing  matrix  equation  for  the  assemblage  of  elements  is  based  on  the  properties 
derived  for  a  single  typical  element.  These  properties,  in  turn,  depend  on  assuming  a  math¬ 
ematical  form  for  the  primary  unknown  of  the  problem  and  then  satisfying  the  variational 
principle.  For  the  elasticity  problem,  the  unknowns  are  the  displacements,  while  for  the  perfect 
incompressible  fluid,  either  the  velocity  potential  or  the  stream  function  may  be  used. 

Of  great  interest  is  the  fact  that  structural  and  nonstructural  elements  may  often  be 
identical  in  shape  and,  further,  be  represented  by  similar  mathematical  expressions.  For 
example,  the  triangle  in  elasticity  will  have  the  same  form  for  displacements  u  and  v,  as  the 
similar  fluid  flow  element  has  for  velocity  potential  4>  or  stream  function  ^ .  In  addition,  the 
numerical  solution  procedure  will  be  quite  similar  for  each  case.  By  way  of  illustration,  the 
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Figure  1.  Saint  Venant  Flow  Problem 


problems  taken  up  in  this  paper  were  solved  using  the  new  ASTRA  structural  program 
developed  at  The  Boeing  Company.  This  versatility,  provided  by  finite  elements  and  associated 
computer  programs,  is  obviously  of  great  practical  importance. 

Finally,  it  should  be  pointed  out  that  the  major  difference  between  the  elasticity  and  fluid 
flow  problems  lies  in  the  boundary  conditions  to  be  satisfied.  This  question  is  therefore 
taken  up  in  the  paper. 
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Figure  2.  Finite  Element  Idealization  of  Cylinder  Flow  Problem 
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SECTION  II 
BASIC  THEORY 


Let  u,  v  represent  the  velocity  components  in  a  two  dimensional  flow.  Then  in  terms  of 
velocity  potential  (x,y) 


U  s 


V  = 


d  4> 
dy 


(I) 


while  in  terms  of  the  stream  function  <//  (x,y) 

d  \j/ 
dy 


v  = 


d\j/ 

dx 


Both  4>  and  ^  must  satisfy  LaPlace’s  equation  or, 


a2  (  )  ,  a2(  ) 

™  T  _ 


=  0 


(2) 


(3) 


dx  dy 

As  shown  in  variational  calculus,  LaPlace’s  equation  may  alternatively  be  expressed  in 
terms  of  a  minimum  principle  (Reference  3).  The  principle  is  as  follows:  that  function  which 
satisfies  LaPlace’s  equation  in  a  region  R  and  which  takes  on  prescribed  values  along  the 
boundary  of  R,  will  have  a  mean  square  gradient  which  is  a  minimum  over  region  R.  Or, 


(4) 


The  equation  also  holds  for  velocity  potential  <f>  .  Note  that  this  is  the  Dirichlet  problem  in 
which  \)/  ,  rather  than  derivatives  of  ^  ,  are  prescribed  on  the  boundary. 


On  solid  boundaries,  values  for  ^  can  be  specified.  At  other  points,  as  the  left  hand 
cross-section  of  Figures  1  or  5,  the  velocity  is  known  and  is  given  in  terms  of  a  derivative 
of  'I'  .  Further  attention  must  be  given  to  such  boundaries.  In  the  case  of  ,  it  is  not  known  on 
solid  walls  or  other  bodies  in  the  flow;  however,  the  flow  d^/dn  normal  to  such  surfaces 
must  be  zero.  In  general  then,  it  will  be  more  difficult  to  treat  the  problem  in  terms  of  <fi  . 
These  points  are  taken  up  in  detail  when  the  examples  of  Figures  1  and  5  are  discussed. 


That  function  ^  (x,y)  which  satisfies  Equation  4,  and  the  boundary  conditions  will  be  the 
solution  to  the  problem.  Such  a  function  must  have  continuous  second  partial  derivatives 
within  each  subregion  of  R  and  must  be  continuous  on  all  element  boundaries.  The  purpose 
will  now  be  to  find  an  approximate  solution  in  terms  of  an  assumed  function  \f/  (x,y). 
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The  assumed  function  is  chosen  with  a  single  element  in  mind.  Such  an  arbitrary  triangle 
is  shown  in  Figure  3.  In  selecting  the  function  f  (x,y),  the  following  points  are  kept  in  mind: 
(1)  the  function  must  be  given  in  terms  of  x  and  y  and  the  values  which  it  has  at  the  nodes 
(1,  2,  3  Figure  3);  (2)  continuity  must  exist  in  f  (x,y)  at  interfaces  of  adjoining  elements; 
(3)  velocities  (Equation  2)  must  include  constant  terms.  These  conditions  will  not  be  dis¬ 
cussed  here.  They  are  necessary  in  order  to  guarantee  boundedness  and  convergence  of  the 
finite  element  solution.  (References  4  and  5). 

Based  on  the  assumed  form  for  \f/  (x,y),  Equation  4  will  be  satisfied  approximately  by 
minimizing  f  with  respect  to  its  nodal  values  j  ^  j.  Or, 

dl  .  o 

»{*, } 

This  will  lead  directly  to  the  fluid  element  stiffness  equation,  which  takes  the  form, 

Mty  H°}  l6> 

where  £  k  j  is  the  fluid  element  stiffness  matrix.  It  will  be  symmetric  and  singular  for 
structural  elements. 


This  equation  for  the  element  may  now  be  treated  precisely  for  the  structural  case  in 
forming  the  total  assemblage  stiffness  equation;  that,  by  a  simple  superposition  procedure 
(Reference  1  or  2).  Let  this  equation  be  given  by. 


["]{*}■{•} 


If  are  unknown,  while  are  prescribed  (the  Dirichlet  problem),  Equation  7  may  be 

3.  D 

written  as, 


Then, 

Koo  *<,  +  “ob  *b  =  0 

*ob  +  “bb  =  0 

Since  is  known,  Equation 9a  may  be  used  to  find  or, 

%  «-  <Kaa>_l  Kab 


(8) 


(9a  ) 

(9b) 


(10) 
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Figure  3.  The  Triangular  Finite  Element 

Note  that  whereas  K  of  Equation  7  is  singular,  the  submatrix  Kaa  of  Equation  9a  will  be 
nonsingular.  Equation  9b  cannot  be  used  in  finding  \fr  since  K T,  need  not  be  a  square 

3  3.0 

matrix  and  hence,  cannot  be  inverted. 

The  derivation  of  the  element  stiffness  matrix,  Equation  6,  is  given  in  the  Appendix. 
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SECTION  m 

SAINT  VENANT’S  PRINCIPLE  IN  FLUID  FLOW 


Saint  Venant’s  principle,  although  widely  known  in  solid  mechanics,  seems  to  be  rather 
infrequently  referred  to  in  fluid  mechanics.  Nevertheless,  this  principle  is  investigated  in 
Reference  6  for  incompressible  potential  flows.  It  is  shown  that  if  the  net  mass  flow  at  a  given 
station  in  a  channel  is  zero,  the  velocity  profile  will  approach  zero  a  short  distance  down¬ 
stream.  The  analogy  with  the  structural  case  is  self-evident. 

Figure  1  shows  a  typical  flow  to  which  St.  Venant’s  principle  must  apply.  The  net  flow 
entering  the  channel  is  zero;  hence,  the  velocities  should  quickly  approach  zero  as  the  fluid 
moves  downstream.  The  flow  may  be  taken  as  two  dimensional. 

Boundary  conditions  are  as  follows: 

\f/  -  0  on  y  =  0  ,  I 

o  =  u  on  x  =  0  (t2) 

o 

Due  to  symmetry,  only  the  upper  half  of  the  channel  need  be  considered. 

The  mixed  boundary  value  problem  represented  by  Equations  11  and  12  can  be  easily 
replaced  by  a  Dirichlet  type  problem.  To  do  this,  it  is  merely  necessary  to  use  the  definition 
of  the  stream  function  and  thereby  replace  u  on  x  =  0  with  i' .  With  Figure  1  in  mind,  it 
follows  that, 

%  --  u0y2,  where  \x  z  \ 

=  i 

Similarly,  it  is  found  that  =  0,  =  -1,  and  ^  =  0.  The  result  at  node  5  agrees  with  the 

statement  made  in  Equation  11.  At  a  downstream  boundary,  say  x  =  1,  f  is  unknown  and 
hence,  will  be  found  as  part  of  the  solution. 

The  finite  element  representation  is  indicated  in  Figure  1.  In  terms  of  the  triangular 
elements  shown,  the  solution  procedure  follows  the  steps  illustrated  by  Equations  8,  9,  and  10. 
Results  were  obtained  for  ^  at  the  nodes  and  for  velocity,  u,  for  each  element.  Due  to  the 
linear  form  for  ^  ,  and  in  view  of  Equation  2,  it  is  seen  that  u  can  be  also  calculated  from 
nodal  values  of  the  stream  function. 
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Finite  element  results  agree  with  Reference  6  and  show  that  the  St.  Venant  effect  does 
indeed  occur.  Decay  in  the  velocity  profile  depends  only  on  x  and  was  discovered  to  be  inde¬ 
pendent  of  y.  Therefore,  results  can  be  shown  for  the  strip  bounded  by  horizontal  lines 
starting  at  nodes  1  and  2.  Figure  4  shows  the  velocity  decay  for  this  strip.  The  curve  applies 
equally  well  to  the  other  similar  strips. 

It  is  interesting  to  note  that  the  St.  Venant  effect  occurs  for  the  inviscid  fluid.  Hence,  it 
reflects  the  intrinsic  properties  of  governing  fluid  equations.  These  are  seen  to  be  reflected 
also  by  the  finite  element  representation  of  the  problem. 


Figure  4.  Velocity  Decay  for  Saint  Venant  Flow 
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SECTION  IV 

FLOW  AROUND  CYLINDER  BETWEEN  PARALLEL  WALLS 


A  simple,  but  nevertheless  interesting  boundary  value  problem  is  shown  in  Figure  5. 
Again  the  flow  may  be  taken  as  two  dimensional;  furthermore,  it  is  assumed  that  the  flow  is 
uniform  at  the  upstream  and  downstream  ends  taken  as  shown  in  Figure  5. 


An  approximation  to  this  problem  is  available  in  fluid  theory  (Reference  7).  This  is 
obtained  by  considering  a  distribution  of  doublets  on  a  vertical  line  in  the  presence  of  a 
horizontal  rectilinear  flow.  Unfortunately,  this  scheme  yields  a  circular  cylinder  only  when 
the  strength  of  the  doublets  becomes  small.  As  a  result,  the  theoretical  solution  obtained 
in  this  manner  does  not  fit  the  ratio  of  cylinder  diameter  to  wall  spacing  shown  in  Figure  5. 
The  author  was  unable  to  locate  a  theoretical  solution  valid  for  any  ratio  of  cylinder  diameter 
to  wall  spacing.  It  therefore  may  be  appropriate  to  point  out  at  this  time  that  the  finite  element 
method  is  in  no  sense  restricted  to  parallel  walls  or  a  circular  cylinder. 


Figure  5.  Flow  Around  Cylinder  Between  Parallel  Walls 


For  the  geometry  shown  in  Figure  5,  the  solution  based  on  the  combination  of  doublets 
and  flow  u^  gives 

uc=2.33u0  U(js|.6u0 

These  will  be  commented  on  later. 
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Boundary  conditions  depend  on  the  region  which  will  be  considered  in  the  finite  element 
analysis.  The  smallest,  and  therefore  most  economical,  region  which  can  be  used  is  a-b-c-d-e, 
as  shown  in  Figure  5.  Due  to  symmetry,  it  is  again  permissible  to  choose, 

^  -  0  along  a-b  and  b-c  (14) 

Then  again,  using  the  definition  of  the  stream  function, 

\f/  =  2  on  d  -  e  (l5) 

Along  a-e ,  values  of  \fr  at  nodal  points  can  be  specified  in  the  same  manner  as  used  for  the 
previous  problem.  Finally  on  c-d,  the  vertical  component  of  velocity  must  vanish.  This 
condition  requires  detailed  attention,  and  this  will  be  taken  up  subsequently.  In  the  meantime, 
it  is  pointed  out  that  the  difficulty  is  by-passed  if  the  total  upper  half  of  the  flow  field  is 
considered.  Boundary  conditions  on  a'e 7  are  then  the  same  as  those  on  ae  and  the  formulation 
becomes  once  again  a  Dirichlet  type  problem.  This  is  the  problem  which  was  solved  by 
finite  elements. 

The  idealization  of  the  flow  field  into  a  finite  element  representation  is  shown  in  Figure  2. 
Reflection  of  this  pattern  in  the  computer  program  established  the  right  hand  side  of  the 
flow  field.  Of  interest  is  the  fact  that  this  same  idealization  had  been  previously  used  for 
analyzing  the  corresponding  flat  sheet  with  a  circular  cut-out.  Therefore,  considerable  input 
data  was  already  available  for  undertaking  the  fluid  flow  problem. 

Again  the  computer  furnished  unknown  values  of  ^  and  the  flow  velocity  in  each  triangular 
element.  As  a  result,  a  complete  picture  of  the  flow  field  is  therefore  made  available.  Of 
maximum  interest  is  the  velocity  distribution  along  c-d,  Figure  5.  This  was  calculated 
directly  from  the  nodal  values  of  V' .  These  values  are  listed  in  Table  1.  In  addition,  the 
y-coordinates  of  nodes  lying  on  c-d  are  listed  in  this  same  table. 


TABLE  1 


Node 

— 

y- Coordinate 

* 

67 

1.00 

0 

68 

1.16 

0.3419 

69 

1.33 

0.6880 

70 

1.54 

1.1072 

71 

1.76 

1.5364 

72 

2.00 

2.0000 
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For  the  point  midway  between  nodes  67  and  68,  the  velocity  is  calculated  as, 

\L  -d, 

r«a  r67 

u  : -  =2.20 

87-68  y68  -  *67 

Similarly  for  mid-points  between  the  other  nodes  on  c-d.  These  velocities  are  plotted  in 
Figure  6.  Extrapolation  then  gives  velocity  at  the  cylinder  and  on  the  wall  as  follows: 

uc  =  2 . 20  u0  ud  =  1. 90  u0  (16) 

These  results  can  now  be  compared  with  those  previously  given  in  Equation  13. 

As  seen,  results  compare  well  at  point  c  on  the  cylinder.  Considerable  discrepancy, 
however,  exists  at  point  d  on  the  wall.  The  finite  element  result  for  u^  is  undoubtedly  close 
to  the  true  value  since  from  Figure  6,  it  is  seen  that  it  fits  the  velocity  profile  which  will 
provide  for  continuity  of  the  flow  across  line  c-d  (average  velocity  must  be  2  at  this  section). 


Figure  6.  Velocity  at  Vertical  Centerline  (c-d,  Figure  4) 
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SECTION  V 

BOUNDARY  CONDITIONS  OF  THE  NEUMANN  TYPE 

As  seen  in  the  preceding  problem,  it  may  be  useful  to  prescribe  boundary  conditions  which 
are  given  as  derivatives  of  p .  Generally,  these  would  be  prescribed  on  some  boundaries, 
while  p  is  prescribed  on  others.  This  mixed  type  boundary  value  problem  causes  no  diffi¬ 
culty  in  itself,  provided  a  scheme  can  be  found  for  specifying  boundary  conditions  associated 
with  first  derivatives  of  p .  Furthermore,  the  procedure  developed  for  handling  such  cases 
will  then  be  directly  applicable  to  flow  problems  formulated  in  terms  of  the  velocity  potential 
9.  In  this  latter  case,  however,  since  all  boundary  conditions  would  be  written  in  terms  of 
d^/dn,  a  further  step  is  required  in  order  to  remove  the  singularity  present  in  the  stiffness 
matrix.  As  is  well  known,  the  solution  to  the  Neumann  type  problem  gives  the  potential  function 
to  the  extent  of  an  arbitrary  constant.  Therefore,  in  the  flow  problem  given  in  terms  of  4>  ,  it 
is  permissible  to  arbitrarily  select  4>  at  one  node  (say  p  =  0  at  a  selected  node).  Imposing 
this  condition  on  the  stiffness  equation  for  the  assemblage  will  then  remove  the  singularity 
in  [k  j  .  At  this  point,  the  boundary  conditions  in  terms  of  d  P/d  n  must  be  specified.  Hence, 
for  the  fluid  problem,  whether  formulated  in  terms  of  p  or  p  ,  the  Neumann  type  boundary 
condition  needs  to  be  considered.  Similar  situations  do  arise  in  some  structural  problems 
but  they  have  only  received  scant  attention. 

Returning  to  the  cylinder  problem,  consider  again  the  situation  along  boundary  c-d. 
Since  v  equals  zero  at  this  line,  Equation  2  requires  that, 

dp 

v  1  "  ~dT  =  °  0n  C_d  (17) 

Two  choices  are  available  for  imposing  this  condition  on  the  problem.  First,  the  element 
derivation  which  is  given  in  the  Appendix  can  be  carried  out  with  the  restriction  that  dp/ dx 
be  equal  to  zero.  This  will  give  an  alternate  form  for  the  element  stiffness  matrix  which  can 
then  be  used  for  elements  having  nodes  on  c-d.  The  disadvantage  is,  of  course,  that  such 
elements  must  then  be  incorporated  into  the  computer  program  and  used  as  need  arises.  In 
addition,  versatility  is  lost  in  that  a  condition  such  as  dp/dn  =  0  on  the  cylinder  is  not 
easily  imposed  in  this  manner.  Nevertheless,  for  certain  special  situations  it  may  be  useful 
to  have  such  an  element  in  the  computer  program. 

The  second  possibility  is  to  consider  the  boundary  conditions  on  the  velocity  as  constraints 
on  the  set  of  equations  expressed  as  the  gross  assemblage  stiffness  equation.  Imposing  these 
constraints  is  then  equivalent  to  a  coordinate  transformation  in  which  the  original  set  of  nodal 
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values  for  ^  are  replaced  by  a  new  set  of  generalized  values  *.  In  carrying  out  this 
transformation,  it  must  be  kept  in  mind  that  (or  )  is  a  scalar  and  not  a  vector. 

For  simplicity,  the  element  will  be  discussed  rather  than  the  assemblage.  The  treatment 
is  essentially  the  same  in  each  instance  since  the  element  stiffness  equation  is  similar  in 
form  to  the  corresponding  equation  for  the  total  assemblage.  Also,  the  development  will  be 
carried  out  in  terms  of  ^  although  it  can  just  as  easily  be  written  in  terms  of  4> . 


From  Figure  7  it  is  seen  that, 


d  <p 

=  d  n 


=  Xu  +  fxv 


and  putting  this  equal  to  zero. 


d  if/  ^ 

v»  =  x_dT^a^ 


(18) 


From  the  development  of  f  (x,y)  as  given  in  the  Appendix,  it  can  be  shown  that, 

d\J/  I 


dx  2A 


<*23  +  *31  ^2  +  *12  ^ 


and, 


d\f/  I 

a7=  Ta  1X32  +  x'3  ^  +  Xzt  ^ 


Substituting  into  Equation  18 

(-  Xy23  +  +  (~Xy31  +A*IS  )^2  +(“Xy|2  +/ix2,)^3=0 


Solving  for  in  terms  of  ^  and  'i'g,  enables  the  following  equation  to  be  written, 


i 

1 

0 

% 

^  y  23  ~/xx32 

^*12 

-  H-  X2I 

(19) 

~Xy3(  +  /ax  ,3 

-^31 

*f-/iXl3 

ll/ 

L  r3  _ 

0 

i 

or, 

{*,  }  =  [*]{**}  1201 
Integral  I  of  the  variational  equation  may  be  written  in  terms  of  [k]  and  j  J7.  j  as 
follows  from  the  Appendix.  Or, 

i  {*,}>]{*■  } 


I  = 
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X  =  cos  8 
H-  =  si  n  8 


Figure  7.  Velocities  at  Solid  Boundary 
Substituting  from  Equation  20  gives, 

(21) 

■  i  (+1  } 

Equation  21  provides  the  stiffness  matrix  K*,  subject  to  the  condition  d  <p/d  n  =  0.  It  was 
obtained  by  using  the  transformation  expressed  by  Equations  19  and  20.  Since  these  trans¬ 
formation  equations  can  be  written  for  all  elements  having  nodes  on  the  boundary  where 
d<£/dn  is  to  be  prescribed,  the  treatment  is  applicable  to  the  gross  stiffness  equation.  The 
solution  then  gives  |  from  which  the  original  values  j  ^  j.  can  be  obtained  from  Equation 
20. 


The  scheme  outlined  above  was  not  used  in  solving  the  cylinder  problem  since  it  was 
not  available  in  the  ASTRA  program  at  the  time  the  solution  was  carried  out. 
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SECTION  VI 
DISCUSSION 


Analysis  of  fluid  problems  by  finite  elements  has  barely  begun.  It  is  expected  that  this 
subject  will  receive  intensive  development  in  the  immediate  future.  Of  great  value  for  solving 
problems  will  be  available  structural  elements  and  associated  computer  programs.  These 
will  permit  analysis  of  two  or  three  dimensional  problems.  Extension  to  more  complex 
classes  of  fluid  flow  can  also  be  expected. 

Advantages  of  the  finite  element  method  over  other  numerical  procedures  has  been  well 
established  in  the  existing  literature  (Reference  2). 
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APPENDIX 


The  conditions  as  previously  set  forth  for  choosing  ^  (x,y)  are  met,  in  the  case  of  the 
triangle,  by  taking 


^(x  ,  y)  =  [l  x  y] 
Writing  nodal  values  of  ^  gives, 


[f  (x  ,  y)]  {a} 


(22) 


^2 


I  X 


'2  y2 


I  X 

L 1  **  y3J 


w 


or, 


{+,}  *  [>-]{-> 

Substituting  from  Equation  23  into  22, 

f  (x  ,y)  =  [f  (x  ,  y  )]  [l]  1 


(23) 


(24) 


Equation  24  gives  ^  in  terms  of  its  nodal  values,  Also  note  that  since  $  (x,y)  is  linear, 
it  must  be  identical  for  adjoining  elements  A  and  B  on  common  edge  i-j.  Hence,  continuity  in 
'I'  along  element  interfaces  is  assured.  From  Equation  24, 


(25) 


and  therefore. 


m-  w  w)T[|HT&rw 


A  similar  expression  holds  for 


Since  is  independent  ofx  and  y,  and  similarly  for  •  the  integrand  of  Equation  4 

will  be  constant. 


Therefore, 


(&f  [£])  N'  W 
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where  A  is  the  area  of  the  triangular  element.  Since  I  is  quadratic  in  <//. ,  applying  Equation  5 
gives, 


or. 


where. 


MM  ■  {«} 


‘  3  2  +  1 

c. 

23 

2 

2 

13  X32  +' 

y23 

y  3  1 

X 13 

+ 

y3. 

2.  X32  + 

y  12 

y  23 

X2. 

X  13 

+  y 

(27) 


(SYMMETRIC) 


(28) 


In  this  last  equation,  =  x^ 


x2,  etc. 


The  fluid  element  stiffness  matrix.  Equation  28,  also  appears  as  a  submatrix  in  the 
corresponding  structural  element  stiffness  matrix.  It  is,  therefore,  possible  to  generate  [k  J 
of  Equation  28  from  appropriate  structural  computer  programs. 
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